% chapter02.tex

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %                                                                           %
 %    YABBY documentation                                                    %
 %    Copyright (C) 2007 Vladimir Likic                                      %
 %                                                                           %
 %    The files in this directory provided under the Creative Commons        %
 %    Attribution-NonCommercial-NoDerivs 2.1 Australia license               %
 %    http://creativecommons.org/licenses/by-nc-nd/2.1/au/                   %
 %    See the file license.txt                                               %
 %                                                                           %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\setcounter{section}{0}

\chapter{Tutorial}

The data files used in this tutorial can be found in docs/data/ 
directory. Starting Yabby in this directory should allow one to
execute all examples given in the tutorial.

\section{Working with sequences}

\index{sequences}

In Yabby sequences are represented as sequence objects. A sequence
object may contain one or more sequences. 

One way to create a sequence object is to load sequences from a
file.  Consider the file 'cad3.fas' which contains three sequences
from the Pfam CAD family in the FASTA format:

\begin{verbatim}
>Q53650_STAAU
YVATGIDYLVILILLFSQVKKGQVKHIWIGQYIGTAIVIGASLLVAQGVVNLIPQQWVIG
LLGLLPLYLGVKIWIKGEEDEDESSILSLFSSGKFNQLFLTMIFIVLASSADDFSIYIPY
FTTLSMSEIFIVTIVFLIMVGVLCYVSYRLASFDFISETIEKYERWIVPIVFIGLGIYIL
FENGTSNALISF
>Q97PJ0_STRPN
YISTSIDYLIILIILFAQLSQNKQKWHIYAGQYLGTGLLVGASLVAAYVVNFVPEEWMVG
LLGLIPIYLGIRFAIVGEDAEEEEEEIIERLEQSKANQLFWTVTLLTIASGGDNLGIYIP
YFASLDWSQTLVALLVFVIGIIIFCEISRVLSSIPLIFETIEKYERIIVPLVFILLGLYI
MYENGTIETFLIV
>P95773_STALU
YIAQALDLLVILLMFFARAKTRKEYRDIYIGQYVGSVALIVISLFFAFVLNYVPEKWILG
LLGLIPIYLGIKVAIYGDSDGEERAKKELNEKGLSKLVGTIAIVTIASCGADNIGLFVPY
FVTLSVTNLLITLFVFLILIFFLVFAAQKLANIPEVGEIVEKFGRWIMAVIYIALGLFII
IENDTIQTILGF
\end{verbatim}

\index{seq\_load}

To load this file in the workspace use the command 'seq\_load':

\begin{verbatim}
yabby> seq_load cad3.fas cad3

 Reading the file 'cad3.fas' ..
 3 sequence(s) found.
\end{verbatim}

Most Yabby commands which create a new object in the workspace
conform to the same pattern:

\begin{verbatim}
COMMAND [ options ] ARGS OBJ_NAME
\end{verbatim}

Where '[ options ]' is where flags and options are passed (if
any), ARGS are the arguments to the command, and OBJ\_NAME is
the name of the object to be created. In the example given
above, there were no options to the 'seq\_load' command; there
was one argument, the name of the file; and the object to be
created was named 'cad3'.

It is possible to inspect objects currently available in the
workspace:

\index{what}

\begin{verbatim}
yabby> what

    object(s)      type
  ------------------------------
    cad3           seq           
\end{verbatim}

This listing shows that one object is currently in the workspace,
of the 'seq' type. It is possible to load the same sequences
under a different name: 

\begin{verbatim}
yabby> seq_load cad3.fas cad3_second

 Reading the file 'cad3.fas' ..
 3 sequence(s) found.

yabby> what

    object(s)      type
  ------------------------------
    cad3           seq           
    cad3_second    seq           
\end{verbatim}

The two objects 'cad3.seq' and 'cad3\_second.seq' are identical.

\index{print}

The command 'print' allows one to output the sequence object:

\begin{verbatim}
yabby> print cad3.seq

>Q53650_STAAU [ Q53650_STAAU ]
YVATGIDYLVILILLFSQVKKGQVKHIWIGQYIGTAIVIGASLLVAQGVVNLIPQQWVIG
LLGLLPLYLGVKIWIKGEEDEDESSILSLFSSGKFNQLFLTMIFIVLASSADDFSIYIPY
FTTLSMSEIFIVTIVFLIMVGVLCYVSYRLASFDFISETIEKYERWIVPIVFIGLGIYIL
FENGTSNALISF
>Q97PJ0_STRPN [ Q97PJ0_STRPN ]
YISTSIDYLIILIILFAQLSQNKQKWHIYAGQYLGTGLLVGASLVAAYVVNFVPEEWMVG
LLGLIPIYLGIRFAIVGEDAEEEEEEIIERLEQSKANQLFWTVTLLTIASGGDNLGIYIP
YFASLDWSQTLVALLVFVIGIIIFCEISRVLSSIPLIFETIEKYERIIVPLVFILLGLYI
MYENGTIETFLIV
>P95773_STALU [ P95773_STALU ]
YIAQALDLLVILLMFFARAKTRKEYRDIYIGQYVGSVALIVISLFFAFVLNYVPEKWILG
LLGLIPIYLGIKVAIYGDSDGEERAKKELNEKGLSKLVGTIAIVTIASCGADNIGLFVPY
FVTLSVTNLLITLFVFLILIFFLVFAAQKLANIPEVGEIVEKFGRWIMAVIYIALGLFII
IENDTIQTILGF
\end{verbatim}

Upon reading the sequences Yabby has taken the first string
in the FASTA comment to be the sequence ID, and the full
comment is re-inserted within the square brackets. In this
case the only comment was the ID string, and this is merely
repeated within the square brackets. 

The command 'print' can send sequences to a file, instead
of printing them in the terminal window:

\begin{verbatim}
yabby> print -f tmp.fasta cad3.seq

 'cad3.seq' written to the file 'tmp.fasta'
\end{verbatim}

Currently there are several useful options of the 'print'
command. The option '-l' causes protein sequences to be
printed in a three-letter format:

\begin{verbatim}
yabby> print -l cad3.seq

>Q53650_STAAU
TYR VAL ALA THR GLY ILE ASP TYR LEU VAL ILE LEU 
ILE LEU LEU PHE SER GLN VAL LYS LYS GLY GLN VAL 
....further output deleted....
\end{verbatim}

The option '-t N' truncates all sequences at N residues: 

\begin{verbatim}
yabby> print -t 10 cad3.seq

>Q53650_STAAU [ Q53650_STAAU ]
YVATGIDYLV
>Q97PJ0_STRPN [ Q97PJ0_STRPN ]
YISTSIDYLI
>P95773_STALU [ P95773_STALU ]
YIAQALDLLV
\end{verbatim}

\index{seq\_info}

The command 'seq\_info' prints additional information about
sequence objects:

\begin{verbatim}
yabby> seq_info cad3

 'cad3' contains 3 sequence(s)
   min number of residues: 192 (sequence 'Q53650_STAAU')
   max number of residues: 193 (sequence 'Q97PJ0_STRPN')
\end{verbatim}

The option '-l' causes the number of residues to be printed
for each sequence: 

\begin{verbatim}
yabby> seq_info -l cad3

 'cad3' contains 3 sequence(s)
  1 -> Q53650_STAAU, 192 residues
  2 -> Q97PJ0_STRPN, 193 residues
  3 -> P95773_STALU, 192 residues
\end{verbatim}

It is often required to select one or more sequences from the sequence
object. The command 'seq\_pick' allows one to select sequences from
their order number or sequence ID. For example:

\index{seq\_pick}

\begin{verbatim}
yabby> seq_pick -n 2 cad3 s2       

 Fetching the sequence 2 ('Q97PJ0_STRPN')
 Saving the extracted sequence as 's2'
\end{verbatim}

The above command has picked the sequence number 2 from the 'cad3'
object, and saved this sequence under the name 's2': 

\begin{verbatim}
yabby> what

    object(s)      type
  ------------------------------
    cad3           seq           
    s2             seq           

yabby> print s2.seq

>Q97PJ0_STRPN [ Q97PJ0_STRPN ]
YISTSIDYLIILIILFAQLSQNKQKWHIYAGQYLGTGLLVGASLVAAYVVNFVPEEWMVG
LLGLIPIYLGIRFAIVGEDAEEEEEEIIERLEQSKANQLFWTVTLLTIASGGDNLGIYIP
YFASLDWSQTLVALLVFVIGIIIFCEISRVLSSIPLIFETIEKYERIIVPLVFILLGLYI
MYENGTIETFLIV
\end{verbatim}

Another useful option is to select a sequence by its ID string:

\begin{verbatim}
yabby> seq_pick -q Q53650_STAAU cad3 s1

 Fetching the sequence 'Q53650_STAAU'
 Saving the extracted sequence as 's1'

yabby> what

    object(s)      type
  ------------------------------
    cad3           seq           
    s1             seq           
    s2             seq           
\end{verbatim}

Suppose that we wanted to take the sequence Q53650\_STAAU, and
extract residues 21-40. The above command would take care of the
first part, while the command 'seq\_strip' could be used to
select a residue range:

\index{seq\_strip}

\begin{verbatim}
yabby> seq_strip 21:40 s1 s1_portion

 's1' contains 1 sequence(s)
 stripping 'Q53650_STAAU'

yabby> print s1_portion.seq

>Q53650_STAAU_21:40 [ Q53650_STAAU ]
KGQVKHIWIGQYIGTAIVIG
\end{verbatim}

The command 'seq\_pattern' allows one to search for a pattern
in a sequence. For example, to search for a pattern 'IDY' use: 

\index{seq\_pattern}

\begin{verbatim}
yabby> seq_pattern IDY cad3

 'IDY' matches in 'Q53650_STAAU'
 'IDY' matches in 'Q97PJ0_STRPN'
 3 sequences examined, 2 match(es) found
\end{verbatim} 

The option '-s NAME' allows one to save the matching sequences
under a the new name:

\begin{verbatim}
yabby> seq_pattern -s IDY_matches IDY cad3

 'IDY' matches in 'Q53650_STAAU'
 'IDY' matches in 'Q97PJ0_STRPN'
 3 sequences examined, 2 match(es) found
 Saving matches as 'IDY_matches'
 [ seq_pattern: 'IDY_matches.seq' exists, overwritten ]
\end{verbatim}

The option '-c' allows one to search the sequence comment for a
pattern, rather than the sequence residues:

\begin{verbatim}
yabby> seq_pattern -c STAA cad3

 'STAA' matches in 'Q53650_STAAU'
 3 sequences examined, 1 match(es) found
\end{verbatim}

The command 'seq\_op' allows one to combine two sequence objects
based on their IDs. This command has the following form:

\index{seq\_op}

\begin{verbatim}
seq_op [ options ] SEQ1_OBJ SEQ2_OBJ RES_OBJ
\end{verbatim}

where SEQ1\_OBJ SEQ2\_OBJ are two sequence objects to be combined 
and RES\_OBJ is the name under which the result will be stored.
This command can calculate the union, intersection, symmetric
difference.

For example, consider the sets of sequences stored in files
'cad.fas' and 'cad3.fas', where the first file contains six 
CAD sequences, and cad3.fas contains three out of six sequences
present in 'cad.fas'.

\begin{verbatim}
yabby> seq_load cad.fas cad

 Reading the file 'cad.fas' ..
 6 sequence(s) found.

yabby> seq_load cad3.fas cad3

 Reading the file 'cad3.fas' ..
 3 sequence(s) found.

yabby> seq_op -i cad cad3 cad_i

 Found 6 sequence(s) in 'cad'
 Found 3 sequence(s) in 'cad3'
  INTERSECTION contains 3 sequence(s)
  DIFFERENCE contains 3 sequence(s)
  UNION contains 6 sequence(s)
 [ Saving INTERSECTION as 'cad_i' ]
 [ seq_op: 'cad_i.seq' exists, overwritten ]
\end{verbatim}

The last command has saved the intersection of objects 'cad.seq' and
'cad3.seq' as a new sequence object, named 'cad\_i'. Since 'cad3.seq'
is simply the subset of 'cad.seq' this intersection is identical to
'cad3.seq'. Options '-u' and '-d' save the union and difference.
It should be noted that the command 'seq\_op' compares only
sequence IDs, {\em it does not compare the sequences themselves}.
Furthermore, unpredictable results may occur if the sequence IDs
are not unique within each set of sequences. 

\index{seq\_comment}

The command 'seq\_comment' modifies the comment of each sequence
to append the sequence number to its ID. To modify the IDs of the
sequences in 'cad3':

\begin{verbatim}
yabby> seq_comment cad3

 Comments modified in 3 sequence(s).

yabby> print cad3.seq

>1-Q53650_STAAU [ Q53650_STAAU ]
YVATGIDYLVILILLFSQVKKGQVKHIWIGQYIGTAIVIGASLLVAQGVVNLIPQQWVIG
LLGLLPLYLGVKIWIKGEEDEDESSILSLFSSGKFNQLFLTMIFIVLASSADDFSIYIPY
FTTLSMSEIFIVTIVFLIMVGVLCYVSYRLASFDFISETIEKYERWIVPIVFIGLGIYIL
FENGTSNALISF
>2-Q97PJ0_STRPN [ Q97PJ0_STRPN ]
YISTSIDYLIILIILFAQLSQNKQKWHIYAGQYLGTGLLVGASLVAAYVVNFVPEEWMVG
LLGLIPIYLGIRFAIVGEDAEEEEEEIIERLEQSKANQLFWTVTLLTIASGGDNLGIYIP
YFASLDWSQTLVALLVFVIGIIIFCEISRVLSSIPLIFETIEKYERIIVPLVFILLGLYI
MYENGTIETFLIV
>3-P95773_STALU [ P95773_STALU ]
YIAQALDLLVILLMFFARAKTRKEYRDIYIGQYVGSVALIVISLFFAFVLNYVPEKWILG
LLGLIPIYLGIKVAIYGDSDGEERAKKELNEKGLSKLVGTIAIVTIASCGADNIGLFVPY
FVTLSVTNLLITLFVFLILIFFLVFAAQKLANIPEVGEIVEKFGRWIMAVIYIALGLFII
IENDTIQTILGF
\end{verbatim}

Since this command will change the sequence IDs, the output
of the command 'seq\_op' would be different. Specifically,
the intersection of the two sets of sequences would contain
no sequences:

\begin{verbatim}
yabby> seq_op -i cad cad3 cad_u

 Found 6 sequence(s) in 'cad'
 Found 3 sequence(s) in 'cad3'
  INTERSECTION contains 0 sequence(s)
  DIFFERENCE contains 9 sequence(s)
  UNION contains 9 sequence(s)
 [ No sequences to save ]

yabby> seq_op -u cad cad3 cad_u

 Found 6 sequence(s) in 'cad'
 Found 3 sequence(s) in 'cad3'
  INTERSECTION contains 0 sequence(s)
  DIFFERENCE contains 9 sequence(s)
  UNION contains 9 sequence(s)
 [ Saving UNION as 'cad_u' ]
\end{verbatim}

\index{seq\_unique}

The command 'seq\_unique' finds sequences present in one set
not present in the other, by {\em comparing sequence strings}.

\begin{verbatim}
yabby> seq_unique cad cad3 cadu

 3 unique sequences found.
 Saving sequences as 'cadu'
\end{verbatim}

Note that in this case the order of sequence objects is important.
This command find the sequences present in the first sequence
object and not present in the second, i.e.:

\begin{verbatim}
yabby> seq_unique cad3 cad tmp

 No unique sequences found.
\end{verbatim}

\index{seq\_compl}

The command 'seq\_compl' calculates the reverse complement of
a DNA sequence:

\begin{verbatim}
yabby> seq_load dna.fas dna

 Reading the file 'dna.fas' ..
 1 sequence(s) found.

yabby> print dna.seq

>chr01 [ chr01 ]
taaccctaaccctaaccctgaccctaaccctaaccctaaccctaaccctaaccagtacac
gcgtacacgtacaagcacccgtacccccagtatacttggacacccgtactcagttatcct
ttttattagtgtacccgcctcttgcacgcatgccacagttcttcagcagaagaacacgca
caatgctctttgataaacgtgcggacatgaaaaaaagggaaaaacgcagctacgtgtgct
gtcgttggtttcacagcgtcaagccgcgtcggtgtaccaaagaggaggtgacccatcgag

yabby> seq_compl dna dna_c

 'dna' contains 1 sequence(s)
 Working on 'chr01'

yabby> print dna_c.seq

>chr01 [ chr01 ]
attgggattgggattgggactgggattgggattgggattgggattgggattggtcatgtg
cgcatgtgcatgttcgtgggcatgggggtcatatgaacctgtgggcatgagtcaatagga
aaaataatcacatgggcggagaacgtgcgtacggtgtcaagaagtcgtcttcttgtgcgt
gttacgagaaactatttgcacgcctgtactttttttccctttttgcgtcgatgcacacga
cagcaaccaaagtgtcgcagttcggcgcagccacatggtttctcctccactgggtagctc
\end{verbatim}

The option '-r' calculates the reverse complement: 

\begin{verbatim}
yabby> seq_compl -r dna dna_cr

 'dna' contains 1 sequence(s)
 Working on 'chr01'

yabby> print dna_cr.seq

>chr01 [ chr01 ]
ctcgatgggtcacctcctctttggtacaccgacgcggcttgacgctgtgaaaccaacgac
agcacacgtagctgcgtttttccctttttttcatgtccgcacgtttatcaaagagcattg
tgcgtgttcttctgctgaagaactgtggcatgcgtgcaagaggcgggtacactaataaaa
aggataactgagtacgggtgtccaagtatactgggggtacgggtgcttgtacgtgtacgc
gtgtactggttagggttagggttagggttagggttagggtcagggttagggttagggtta
\end{verbatim}

\index{seq\_genbank}

The command 'seq\_genbank' fetches a sequence from GenBank
by default using the sequence accession number:

\begin{verbatim}
yabby> seq_genbank -a J00522 mig

 Saving 'J00522' as 'mig'

yabby> what

    object(s)      type
  ------------------------------
    mig            seq           
\end{verbatim}

This command requires the module Bio::DB::GenBank and illustrates how
BioPerl \cite{bioperl} can be used within Yabby.

\index{seq\_os}

The command 'seq\_os' inserts additional information into sequence
comments, where additional information is read from an external
file. For example, the file 'sprot\_test.dat.os' contains the
list of sequences IDs and associated organism names, one line
per sequence:

\begin{verbatim}
104K_THEAN Theileria annulata.
104K_THEPA Theileria parva.
108_SOLLC Solanum lycopersicum (Tomato) (Lycopersicon esculentum).
10KD_VIGUN Vigna unguiculata (Cowpea).
110KD_PLAKN Plasmodium knowlesi.
11S2_SESIN Sesamum indicum (Oriental sesame) (Gingelly).
11S3_HELAN Helianthus annuus (Common sunflower).
11SB_CUCMA Cucurbita maxima (Pumpkin) (Winter squash).
128UP_DROME Drosophila melanogaster (Fruit fly).
12AH_CLOS4 Clostridium sp. (strain C 48-50).
12KD_FRAAN Fragaria ananassa (Strawberry).
12KD_MYCSM Mycobacterium smegmatis.
12S1_ARATH Arabidopsis thaliana (Mouse-ear cress).
12S2_ARATH Arabidopsis thaliana (Mouse-ear cress).
12S_PROFR Propionibacterium freudenreichii subsp. shermanii.
13KDA_SCYCA Scyliorhinus canicula (Spotted dogfish) (Spotted catshark).
13KDA_TRISC Triakis scyllium (Leopard shark) (Triakis scyllia).
13S1_FAGES Fagopyrum esculentum (Common buckwheat).
\end{verbatim}

The command 'seq\_os' could be used to insert the organism name
in the sequence object comments by using the file listed above.
To illustrate this, we first load the test set of sequences:

\begin{verbatim}
yabby> seq_load sprot_test.fas test

 Reading the file 'sprot_test.fas' ..
 3 sequence(s) found.

yabby> print test.seq

>104K_THEPA [ 104K_THEPA 104 kDa microneme/rhoptry antigen precursor (p104) ]
MKFLILLFNILCLFPVLAADNHGVGPQGASGVDPITFDINSNQTGPAFLTAVEMAGVKYL
QVQHGSNVNIHRLVEGNVVIWENASTPLYTGAIVTNNDGPYMAYVEVLGDPNLQFFIKSG
DAWVTLSEHEYLAKLQEIRQAVHIESVFSLNMAFQLENNKYEVETHAKNGANMVTFIPRN
GHICKMVYHKNVRIYKATGNDTVTSVVGFFRGLRLLLINVFSIDDNGMMSNRYFQHVDDK
YVPISQKNYETGIVKLKDYKHAYHPVDLDIKDIDYTMFHLADATYHEPCFKIIPNTGFCI
TKLFDGDQVLYESFNPLIHCINEVHIYDRNNGSIICLHLNYSPPSYKAYLVLKDTGWEAT
THPLLEEKIEELQDQRACELDVNFISDKDLYVAALTNADLNYTMVTPRPHRDVIRVSDGS
EVLWYYEGLDNFLVCAWIYVSDGVASLVHLRIKDRIPANNDIYVLKGDLYWTRITKIQFT
QEIKRLVKKSKKKLAPITEEDSDKHDEPPEGPGASGLPPKAPGDKEGSEGHKGPSKGSDS
SKEGKKPGSGKKPGPAREHKPSKIPTLSKKPSGPKDPKHPRDPKEPRKSKSPRTASPTRR
PSPKLPQLSKLPKSTSPRSPPPPTRPSSPERPEGTKIIKTSKPPSPKPPFDPSFKEKFYD
DYSKAASRSKETKTTVVLDESFESILKETLPETPGTPFTTPRPVPPKRPRTPESPFEPPK
DPDSPSTSPSEFFTPPESKRTRFHETPADTPLPDVTAELFKEPDVTAETKSPDEAMKRPR
SPSEYEDTSPGDYPSLPMKRHRLERLRLTTTEMETDPGRMAKDASGKPVKLKRSKSFDDL
TTVELAPEPKASRIVVDDEGTEADDEETHPPEERQKTEVRRRRPPKKPSKSPRPSKPKKP
KKPDSAYIPSILAILVVSLIVGIL
>110KD_PLAKN [ 110KD_PLAKN 110 kDa antigen (PK110) (Fragment) ]
FNSNMLRGSVCEEDVSLMTSIDNMIEEIDFYEKEIYKGSHSGGVIKGMDYDLEDDENDED
EMTEQMVEEVADHITQDMIDEVAHHVLDNITHDMAHMEEIVHGLSGDVTQIKEIVQKVNV
AVEKVKHIVETEETQKTVEPEQIEETQNTVEPEQTEETQKTVEPEQTEETQNTVEPEQIE
ETQKTVEPEQTEEAQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTE
ETQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTEETQNTVEPEPTQETQNTVEP
>12KD_FRAAN [ 12KD_FRAAN Auxin-repressed 12.5 kDa protein. ]
MVLLDKLWDDIVAGPQPERGLGMLRKVPQPLNLKDEGESSKITMPTTPTTPVTPTTPISA
RKDNVWRSVFHPGSNLSSKTMGNQVFDSPQPNSPTVYDWMYSGETRSKHHR
\end{verbatim}

The command 'seq\_os' takes two arguments, the file name with 
sequence IDs and additional information (see above), and the
sequence object name:

\begin{verbatim}
yabby> seq_os sprot_test.dat.os test

 Reading the file 'sprot_test.dat.os' ...
 3 sequences to process.

 All done.
 [ seq_os: 'test.seq' exists, overwritten ]
\end{verbatim}

The original sequence object was overwritten, but the difference
is only in that the comments are changed to contain the information
provided in the file 'sprot\_test.dat.os':

\begin{verbatim}

yabby> print test.seq

>104K_THEPA [ 104K_THEPA 104 kDa microneme/rhoptry antigen
... precursor (p104) Theileria parva. ]
MKFLILLFNILCLFPVLAADNHGVGPQGASGVDPITFDINSNQTGPAFLTAVEMAGVKYL
QVQHGSNVNIHRLVEGNVVIWENASTPLYTGAIVTNNDGPYMAYVEVLGDPNLQFFIKSG
DAWVTLSEHEYLAKLQEIRQAVHIESVFSLNMAFQLENNKYEVETHAKNGANMVTFIPRN
GHICKMVYHKNVRIYKATGNDTVTSVVGFFRGLRLLLINVFSIDDNGMMSNRYFQHVDDK
YVPISQKNYETGIVKLKDYKHAYHPVDLDIKDIDYTMFHLADATYHEPCFKIIPNTGFCI
TKLFDGDQVLYESFNPLIHCINEVHIYDRNNGSIICLHLNYSPPSYKAYLVLKDTGWEAT
THPLLEEKIEELQDQRACELDVNFISDKDLYVAALTNADLNYTMVTPRPHRDVIRVSDGS
EVLWYYEGLDNFLVCAWIYVSDGVASLVHLRIKDRIPANNDIYVLKGDLYWTRITKIQFT
QEIKRLVKKSKKKLAPITEEDSDKHDEPPEGPGASGLPPKAPGDKEGSEGHKGPSKGSDS
SKEGKKPGSGKKPGPAREHKPSKIPTLSKKPSGPKDPKHPRDPKEPRKSKSPRTASPTRR
PSPKLPQLSKLPKSTSPRSPPPPTRPSSPERPEGTKIIKTSKPPSPKPPFDPSFKEKFYD
DYSKAASRSKETKTTVVLDESFESILKETLPETPGTPFTTPRPVPPKRPRTPESPFEPPK
DPDSPSTSPSEFFTPPESKRTRFHETPADTPLPDVTAELFKEPDVTAETKSPDEAMKRPR
SPSEYEDTSPGDYPSLPMKRHRLERLRLTTTEMETDPGRMAKDASGKPVKLKRSKSFDDL
TTVELAPEPKASRIVVDDEGTEADDEETHPPEERQKTEVRRRRPPKKPSKSPRPSKPKKP
KKPDSAYIPSILAILVVSLIVGIL
>110KD_PLAKN [ 110KD_PLAKN 110 kDa antigen (PK110) (Fragment)
... Plasmodium knowlesi. ]
FNSNMLRGSVCEEDVSLMTSIDNMIEEIDFYEKEIYKGSHSGGVIKGMDYDLEDDENDED
EMTEQMVEEVADHITQDMIDEVAHHVLDNITHDMAHMEEIVHGLSGDVTQIKEIVQKVNV
AVEKVKHIVETEETQKTVEPEQIEETQNTVEPEQTEETQKTVEPEQTEETQNTVEPEQIE
ETQKTVEPEQTEEAQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTE
ETQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTEETQNTVEPEPTQETQNTVEP
>12KD_FRAAN [ 12KD_FRAAN Auxin-repressed 12.5 kDa protein.
... Fragaria ananassa (Strawberry). ]
MVLLDKLWDDIVAGPQPERGLGMLRKVPQPLNLKDEGESSKITMPTTPTTPVTPTTPISA
RKDNVWRSVFHPGSNLSSKTMGNQVFDSPQPNSPTVYDWMYSGETRSKHHRim}

\end{verbatim}

\subsection{Swiss-Prot related commands}

Several commands exist for the manipulation of Swiss-Prot data
files.

\index{sprot\_fetch}

The command 'sprot\_fetch' fetches a single sequence from the
Swiss-Prot data file based on the sequence unique ID. For example,
the Swiss-Prot file 'sprot\_test.dat' contains 18 sequences,
including the sequence with the ID '110KD\_PLAKN'. The command
'sprot\_fetch' could be used to extract this sequence from the
database file:

\begin{verbatim}
yabby> sprot_fetch sprot_test.dat 110KD_PLAKN plakn

 Sequence '110KD_PLAKN' found.
 [ Saving as 'plakn.seq' ]

yabby> print plakn.seq

>110KD_PLAKN [ 110 kDa antigen (PK110) (Fragment). Plasmodium knowlesi. ]
FNSNMLRGSVCEEDVSLMTSIDNMIEEIDFYEKEIYKGSHSGGVIKGMDYDLEDDENDED
EMTEQMVEEVADHITQDMIDEVAHHVLDNITHDMAHMEEIVHGLSGDVTQIKEIVQKVNV
AVEKVKHIVETEETQKTVEPEQIEETQNTVEPEQTEETQKTVEPEQTEETQNTVEPEQIE
ETQKTVEPEQTEEAQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTE
ETQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTEETQNTVEPEPTQETQNTVEP
\end{verbatim}

\index{seq\_sprot\_os}

The command 'seq\_sprot\_os' fetches the organism name from the
Swiss-Prot data file (Swiss-Prot field "OS") for a set of sequences
saved in the workspace. This is useful when the organism name is
lost in the sequence manipulation pipeline. For example, when the
Swiss-Prot database is convereted into the FASTA format with the
program "sreformat" from the Sean Eddy's SQUID toolkit, the organism
names are not copied to the FASTA file.

The command 'seq\_sprot\_os' requires two arguments, the name of the
Swiss-Prot database file and the name of an existing sequence object.
To illustrate this, we first load the test set of sequences
('sprot\_test.fas') in the workspace, under the name 'test':

\begin{verbatim}
yabby> seq_load sprot_test.fas test 

 Reading the file 'sprot_test.fas' ..
 3 sequence(s) found.

yabby> print test.seq

>104K_THEPA [ 104K_THEPA 104 kDa microneme/rhoptry antigen precursor (p104) ]
MKFLILLFNILCLFPVLAADNHGVGPQGASGVDPITFDINSNQTGPAFLTAVEMAGVKYL
QVQHGSNVNIHRLVEGNVVIWENASTPLYTGAIVTNNDGPYMAYVEVLGDPNLQFFIKSG
DAWVTLSEHEYLAKLQEIRQAVHIESVFSLNMAFQLENNKYEVETHAKNGANMVTFIPRN
GHICKMVYHKNVRIYKATGNDTVTSVVGFFRGLRLLLINVFSIDDNGMMSNRYFQHVDDK
YVPISQKNYETGIVKLKDYKHAYHPVDLDIKDIDYTMFHLADATYHEPCFKIIPNTGFCI
TKLFDGDQVLYESFNPLIHCINEVHIYDRNNGSIICLHLNYSPPSYKAYLVLKDTGWEAT
THPLLEEKIEELQDQRACELDVNFISDKDLYVAALTNADLNYTMVTPRPHRDVIRVSDGS
EVLWYYEGLDNFLVCAWIYVSDGVASLVHLRIKDRIPANNDIYVLKGDLYWTRITKIQFT
QEIKRLVKKSKKKLAPITEEDSDKHDEPPEGPGASGLPPKAPGDKEGSEGHKGPSKGSDS
SKEGKKPGSGKKPGPAREHKPSKIPTLSKKPSGPKDPKHPRDPKEPRKSKSPRTASPTRR
PSPKLPQLSKLPKSTSPRSPPPPTRPSSPERPEGTKIIKTSKPPSPKPPFDPSFKEKFYD
DYSKAASRSKETKTTVVLDESFESILKETLPETPGTPFTTPRPVPPKRPRTPESPFEPPK
DPDSPSTSPSEFFTPPESKRTRFHETPADTPLPDVTAELFKEPDVTAETKSPDEAMKRPR
SPSEYEDTSPGDYPSLPMKRHRLERLRLTTTEMETDPGRMAKDASGKPVKLKRSKSFDDL
TTVELAPEPKASRIVVDDEGTEADDEETHPPEERQKTEVRRRRPPKKPSKSPRPSKPKKP
KKPDSAYIPSILAILVVSLIVGIL
>110KD_PLAKN [ 110KD_PLAKN 110 kDa antigen (PK110) (Fragment) ]
FNSNMLRGSVCEEDVSLMTSIDNMIEEIDFYEKEIYKGSHSGGVIKGMDYDLEDDENDED
EMTEQMVEEVADHITQDMIDEVAHHVLDNITHDMAHMEEIVHGLSGDVTQIKEIVQKVNV
AVEKVKHIVETEETQKTVEPEQIEETQNTVEPEQTEETQKTVEPEQTEETQNTVEPEQIE
ETQKTVEPEQTEEAQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTE
ETQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTEETQNTVEPEPTQETQNTVEP
>12KD_FRAAN [ 12KD_FRAAN Auxin-repressed 12.5 kDa protein. ]
MVLLDKLWDDIVAGPQPERGLGMLRKVPQPLNLKDEGESSKITMPTTPTTPVTPTTPISA
RKDNVWRSVFHPGSNLSSKTMGNQVFDSPQPNSPTVYDWMYSGETRSKHHR
\end{verbatim}

The Swiss-Prot example file 'sprot\_test.dat' contains 18
sequences which include the tree from the file 'sprot\_test.fas'.
The command 'seq\_sprot\_os' could be used to fetch the OS
field from the  Swiss-Prot file, and insert this in the
sequence comment:

\begin{verbatim}
yabby> seq_sprot_os sprot_test.dat test

 3 sequences to process.
 Printing dot per processed sequence:
 ...
 All done.
 [ seq_sprot_os: 'test.seq' exists, overwritten ]
\end{verbatim}

This manipulation has happened in place, and the original
sequence object was overwritten. the difference is of course
only in the comment field:

\begin{verbatim}
yabby> print test.seq

>104K_THEPA [ 104K_THEPA 104 kDa microneme/rhoptry antigen precursor
... (p104) Theileria parva. ]
MKFLILLFNILCLFPVLAADNHGVGPQGASGVDPITFDINSNQTGPAFLTAVEMAGVKYL
QVQHGSNVNIHRLVEGNVVIWENASTPLYTGAIVTNNDGPYMAYVEVLGDPNLQFFIKSG
DAWVTLSEHEYLAKLQEIRQAVHIESVFSLNMAFQLENNKYEVETHAKNGANMVTFIPRN
GHICKMVYHKNVRIYKATGNDTVTSVVGFFRGLRLLLINVFSIDDNGMMSNRYFQHVDDK
YVPISQKNYETGIVKLKDYKHAYHPVDLDIKDIDYTMFHLADATYHEPCFKIIPNTGFCI
TKLFDGDQVLYESFNPLIHCINEVHIYDRNNGSIICLHLNYSPPSYKAYLVLKDTGWEAT
THPLLEEKIEELQDQRACELDVNFISDKDLYVAALTNADLNYTMVTPRPHRDVIRVSDGS
EVLWYYEGLDNFLVCAWIYVSDGVASLVHLRIKDRIPANNDIYVLKGDLYWTRITKIQFT
QEIKRLVKKSKKKLAPITEEDSDKHDEPPEGPGASGLPPKAPGDKEGSEGHKGPSKGSDS
SKEGKKPGSGKKPGPAREHKPSKIPTLSKKPSGPKDPKHPRDPKEPRKSKSPRTASPTRR
PSPKLPQLSKLPKSTSPRSPPPPTRPSSPERPEGTKIIKTSKPPSPKPPFDPSFKEKFYD
DYSKAASRSKETKTTVVLDESFESILKETLPETPGTPFTTPRPVPPKRPRTPESPFEPPK
DPDSPSTSPSEFFTPPESKRTRFHETPADTPLPDVTAELFKEPDVTAETKSPDEAMKRPR
SPSEYEDTSPGDYPSLPMKRHRLERLRLTTTEMETDPGRMAKDASGKPVKLKRSKSFDDL
TTVELAPEPKASRIVVDDEGTEADDEETHPPEERQKTEVRRRRPPKKPSKSPRPSKPKKP
KKPDSAYIPSILAILVVSLIVGIL
>110KD_PLAKN [ 110KD_PLAKN 110 kDa antigen (PK110) (Fragment)
... Plasmodium knowlesi. ]
FNSNMLRGSVCEEDVSLMTSIDNMIEEIDFYEKEIYKGSHSGGVIKGMDYDLEDDENDED
EMTEQMVEEVADHITQDMIDEVAHHVLDNITHDMAHMEEIVHGLSGDVTQIKEIVQKVNV
AVEKVKHIVETEETQKTVEPEQIEETQNTVEPEQTEETQKTVEPEQTEETQNTVEPEQIE
ETQKTVEPEQTEEAQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTE
ETQKTVEPEQTEETQKTVEPEQTEETQKTVEPEQTEETQNTVEPEPTQETQNTVEP
>12KD_FRAAN [ 12KD_FRAAN Auxin-repressed 12.5 kDa protein.
... Fragaria ananassa (Strawberry). ]
MVLLDKLWDDIVAGPQPERGLGMLRKVPQPLNLKDEGESSKITMPTTPTTPVTPTTPISA
RKDNVWRSVFHPGSNLSSKTMGNQVFDSPQPNSPTVYDWMYSGETRSKHHR
\end{verbatim}

Note that all sequences from the sequence object must be
present in the Swiss-Prot data file, or this command will
fail (in the example above, all sequences from the 'test.seq'
object must be present in the file 'sprot\_test.dat'). On
the other hand, the Swiss-Prot file can contain additional
sequences. The sequences are searched by the sequence ID. 

The usability of commands 'sprot\_fetch' and 'seq\_sprot\_os'
is limited if the Swiss-Prot file is large and/or the sequence
object is large. For example, the UniProt release 12.4
(23-Oct-2007) contains 5,275,429 sequences, and the combined
Swiss-Prot file is approx 13 Gb in size (swiss-prot + TrEMBL).
The command 'seq\_sprot\_os' cannot cope well with the file
of this size.

\index{sprot\_os}

Consider for example a sequence object Omp85 with 1,000
sequences, which was derived in some way from the UniProt
database (for example, the list of hits from the hidden
Markov model search of the combined UniProt data file).
To associate organism names with the Omp85 sequences, one
would need to search the entire UniProt for each of the
1,000 sequences in the Omp85.  The command 'seq\_sprot\_os'
is not optimized to cope with such a large computational 
task. However, this task can be accomplished in two steps,
by using the commands 'sprot\_os' and 'seq\_os'. First,
the command 'sprot\_os' can be used to extract the organism
names from the UniProt file:

\begin{verbatim}
yabby> sprot_os uniprot.dat uniprot.dat.os

 Processing the database file 'uniprot.dat' ..
 Printing a dot per 10000 processed sequences:
 ............................................................
 ............................................................
 ............................................................
 ............................................................
 ............................................................
 ............................................................
 ............................................................
 ............................................................
 ...............................................

 All done.
\end{verbatim}

In the above example, 'uniprot.dat' is the combined UniProt
data file (swiss-prot + TrEMBL) in the Swiss-Prot format.
The organism names are saved in the file 'uniprot.dat.os'.
This file lists sequence IDs and the first OS field as given
in the Swiss-Prot file, one line per sequence:

\begin{verbatim}
A0AQI4_9ARCH uncultured archaeon.
A0AQI5_9ARCH uncultured archaeon.
A0AQI7_9ARCH uncultured archaeon.
A0AQI8_9ARCH uncultured archaeon.
A0B530_METTP Methanosaeta thermophila (strain DSM 6194 / PT)
A0B531_METTP Methanosaeta thermophila (strain DSM 6194 / PT)
A0B532_METTP Methanosaeta thermophila (strain DSM 6194 / PT)
A0B533_METTP Methanosaeta thermophila (strain DSM 6194 / PT)
A0B534_METTP Methanosaeta thermophila (strain DSM 6194 / PT)
A0B535_METTP Methanosaeta thermophila (strain DSM 6194 / PT)
...
\end{verbatim}

The command 'seq\_os' (explained above) can be used to insert
the organism names in the sequence object, by using the
intermendiate file created by the command 'sprot\_os':

\begin{verbatim}
yabby> seq_os uniprot.dat.os Omp85

 Reading the file 'uniprot.dat.os' ...
 1000 sequences to process.

 All done.
 [ seq_os: 'Omp85.seq' exists, overwritten ]
\end{verbatim}

Similarly as with the command 'sprot\_test.dat', the insertion
of comments was in place, and the sequence object 'Omp85.seq'
was overwritten.  

\section{Housekeeping commands}

\index{help}

Several commands are dedicated to the general manipulation of
Yabby objects. These include 'what' and 'print', and already
seen 'help, 'dump', 'restore', 'delete', 'flush', 'license',
'exit' and 'quit'.

The 'help' command is probably the most useful. By itself it
prints the list of currently available commands: 

\begin{verbatim}
yabby> help

 39 commands ready.

 blast          blast_info     blastg         delete
 dump           emboss_needle  exit           flush
 help           hmm_score2seq  hmm_score_extract license
 mol2seq        mol_crd        mol_load       motif_cmp
 motif_load     motif_meme     pdb_conv       pdb_model
 pfam_fetch     print          quit           restore
 seq2id         seq_comment    seq_compl      seq_genbank
 seq_info       seq_letter     seq_letterc    seq_load
 seq_op         seq_pattern    seq_pick       seq_strip
 seq_unique     sprot_split    what

 For info about a particular command try: 'help COMMAND'
\end{verbatim}

If followed by a command name, it prints the information
about the command:

\begin{verbatim}
yabby> help seq_load


 Loads sequence(s) from the database file.

 Usage:
        seq_load [ options ] DBA_FILE OBJ_NAME

 Where DBA_FILE is the name of the database file. OBJ_NAME is
 the internal Yabby name for the sequence(s).

 Options:

 -a -- append sequences to an already existing sequence
  object OBJ_NAME
 -f -- the file format is FASTA (default)
 -b -- the file format is BLOCKS

 Notes:
 1. Only BLOCKS format as given by MEME output was tested.
\end{verbatim}

\index{delete}
\index{flush}

The command 'delete' deletes an object, and the command 'flush'
deletes all objects from the workspace:

\begin{verbatim}
yabby> what

    object(s)      type
  ------------------------------
    hits           seq           
                   hmm_score     

yabby> delete hits.seq

 [ 'hits.seq' deleted ]

yabby> flush

 [ workspace flushed ]

yabby> what

    object(s)      type
  ------------------------------
   [ none ]
\end{verbatim}

\index{dump}
\index{restore}

The command 'dump' stores the session onto a file:

\begin{verbatim}
yabby> dump mysession

 Yabby session archived as 'mysession.tar.gz'
\end{verbatim}

The session can be restored with the command 'restore':

\begin{verbatim}
yabby> restore mysession

 Yabby session 'mysession' restored
\end{verbatim}

Note that the commands 'dump'/'restore' depend on unix commands
'tar' and 'gzip' being available for execution by Perl.

\section{Working with sequence motifs}

\index{sequence motifs}

Sequence motifs are handled as 'motif' objects.  Motifs could
be loaded from FASTA files, or from BLOCKS files:

\index{motif\_load}

\begin{verbatim}
yabby> motif_load -b m2.blocks m2

 Reading the file 'm2.blocks' ..
 11 sequence(s) found in the motif 'm2'.
\end{verbatim}

Motifs can be printed:

\begin{verbatim}
yabby> what      

    object(s)      type
  ------------------------------
    m2             motif         

yabby> print m2.motif

>Sthermophile [ from BLOCKS file 'm2.blocks' ]
PMAESGVTIRSDSEQYSRHEEQSVSPPSSSS
>Panserina [ from BLOCKS file 'm2.blocks' ]
PMAESGVTIRSDSEQYSSPEELSTSPPSSSS
>Cglobosum [ from BLOCKS file 'm2.blocks' ]
PMAESGITIPSDSEQYSRHGDQSTSPPSSSS
>Ncrassa [ from BLOCKS file 'm2.blocks' ]
PLAESGVTISSDSEQYSAPESASPQSPSSSS
>Oclavigerum [ from BLOCKS file 'm2.blocks' ]
PMAESGVTMHSDSEQYSGAEELSASLESSHS
>Afumigatus [ from BLOCKS file 'm2.blocks' ]
ELYDSGLTVPSDSEIYSANHEVSSPMSASSS
>Gzea [ from BLOCKS file 'm2.blocks' ]
HLAESGVTMHSDIELYSAGDDLSSPPSSNSP
>Aoryzae [ from BLOCKS file 'm2.blocks' ]
ELYESGLTVRSDSENYSANNELSESTSSSPL
>Tlanuginosus [ from BLOCKS file 'm2.blocks' ]
ELSDSALTVPSDSENYSANNEFSSSPSASNS
>Pnodorum [ from BLOCKS file 'm2.blocks' ]
PLTASGLTIPTDSESYSAPADSPSPSPPSSS
>Apullulans [ from BLOCKS file 'm2.blocks' ]
RDIYDSMTMQSDSETYDQPDISSPSSPSSDS
\end{verbatim}

\index{motif\_cmp}

The command 'motif\_cmp' compares two motifs. Two motifs are
identical if they have the same number of sequences, the sequences
have the same ID, and the sequences themselves are identical
as strings:

\begin{verbatim}
yabby> motif_cmp m2 m2_second

 Motifs 'm2' and 'm2_second' contain the same sequence IDs.
 Comparing the sequences...
 Motifs 'm2' and 'm2_second' are identical.
\end{verbatim}

\index{motif\_meme}
\index{MEME}

The command 'motif\_meme' extracts the motifs from the plain
text output of the motif detection program MEME \cite{meme}.
To extract the motif 1 from the file 'meme.out' use:

\begin{verbatim}
yabby> motif_meme meme.out 1 m1 

 Reading MEME output 'meme.out' ..
 Motif 1 saved as 'm1'.

yabby> what

    object(s)      type
  ------------------------------
    m1             motif         
\end{verbatim}

In the internal representation 'motif' type object is identical to
the 'seq' type object.

\section{Working with the HMMER output}

\index{HMMER}
\index{HMMPFAM}
\index{hmm\_score\_extract}
\index{hmm\_score2seq}

HMMER is a powerful software for the profile hidden Markov model
search of sequence databases \cite{hmmer,wwwhmmer}. The 'search 
step' in HMMER search pipeline is performed by the program
HMMPFAM which is a part of the HMMER software suite. Yabby's
commands 'hmm\_score\_extract' and 'hmm\_score2seq' can aid in the
extraction of the search results from the HMMPFAM output files.

As an example, consider the HMMPFAM search performed
on the genome sequence 'LmjFmockup.pep'. The output of this search
was redirected to the file 'hmmpfam.out'. To extract the hits
with the E-value above 0.01, and store these hits under the name
'hits':

\begin{verbatim}
yabby> hmm_score_extract -E 0.01 -s hits hmmpfam.out

 Processing HMMPFAM search output file

  No     Sequence                            E-score
 -----------------------------------------------------
 (   1)  LmjF05.1190@All                     2.70e-03
 (   2)  LmjF05.1190@GlycogenStarch          4.20e-03
 (   3)  LmjF05.0920@All                     7.50e-03
yabby> what

    object(s)      type
  ------------------------------
    hits           hmm_score     
\end{verbatim}

The command 'hmm\_score\_extract' with the '-s hits' flag has created
the object of the type 'hmm\_score' in the workspace. This is
merely a list containing the sequence name, the model the
hit was recorded to, and the E-value:

\begin{verbatim}
LmjF05.1190@All 0.0027
LmjF05.1190@GlycogenStarch 0.0042
LmjF05.0920@All 0.0075
\end{verbatim}

To get to the actual sequences, one needs to go to the
original sequence database on which the HMMPFAM search
was executed, match the sequence IDs, and extract the
sequences. This can be achieved with the command
'hmm\_score2seq'

\begin{verbatim}
yabby> hmm_score2seq LmjFmockup.pep hits.fas hits

 found 3 sequences to extract
 Processing the database 'LmjFmockup.pep'
 Sequences written to 'hits.fas'
\end{verbatim}

This command takes the name of the sequence database on
which the search was executed (LmjFmockup.pep), the name
of the output file where sequence are to be saved (hits.fas),
and the name of the Yabby 'hmm\_score' object where the
hits were extracted.

The output file contains sequence hits in the FASTA format,
which can be loaded and manipulated with sequence commands: 

\begin{verbatim}
yabby> seq_load hits.fas  hits

 Reading the file 'hits.fas' ..
 3 sequence(s) found.

yabby> print hits.seq

>LmjF05.1190@All [ LmjF05.1190@All [ (score: 2.70e-03, model: All)
MSNKKQTWANAHSQLWTSSVARSAGSRTQSEVSSIASTNRDRSLLDDGDGYHQPPSHVCA
QELRHQAYQPNDKTLAQVSDILESYGAIACDADIPAVLLEVVKELERTKRDNNFKDLLLR
EYSDTVQRRFGLYGEESPPSIAEVSARLRDGAPARPVAAPNPAWLEQPLNELWRTVCDGM
QACFEKNADLSAPVLLPEARRTKANVSQLLHTSCDVLSQLAKEYTTAKAAVMKKMERTVA
ASQSFRQLTLSTALSELEAQQELKGSAAEDGGASAMDGSHTVGSALQSLIDVIPASLQIH
.... extra output deleted ....
\end{verbatim}


\section{Running NCBI BLAST from Yabby}

BLAST (Basic Local Alignment Search Tool) refers to a heuristic
algorithm widely used to find regions of similarity between
protein/DNA sequences. BLAST is often used to denote the computer
program, an actual implementation of the BLAST algorithm. Two
such implementations are in wide use: NCBI BLAST and WU-BLAST
(from the Washington University). In addition, the term BLAST
is often used to denote the Web interface for blast search,
such as the one provided by NCBI \cite{wwwblast}.

It is generally less known that NCBI BLAST executables can be
downloaded \cite{ftpblast} and run locally, which has some
distinct advantages. The Yabby command 'blast' allows one to
run NCBI BLAST search locally, through the Yabby interface. 

The command 'blast' depends on NCBI BLAST being installed locally.
The location of the NCBI BLAST installation is set in the
Yabby library 'yabby\_blast.pm':

\begin{verbatim}
$BLASTALL = "/usr/local/bin/blastall";
\end{verbatim}

For NCBI BLAST to work one also needs to set up the file
'.ncbirc' in the user's home directory that will give the
location of blast data files (such as substitution matrices): 

\begin{verbatim}
[NCBI]
   Data=/usr/local/blast/2.2.9/data
\end{verbatim}

The database to be searched need first to be indexed for BLAST
search with the program 'formatdb' which is provided in the
NCBI BLAST package.  For example, to format the database
'LmjFwholegenome.pep' (L.major genome in FASTA fomat),

\begin{verbatim}
$ /usr/local/bin/formatdb -i LmjFwholegenome.pep -p T -o T
\end{verbatim}

This would create several additional files in the same directory,
such as 'LmjFwholegenome.pep.phr', 'LmjFwholegenome.pep.psq' etc.

Yabby can execute BLAST search against multiple sequences. For
example, to BLAST six sequences from the file 'cad.fas':

\begin{verbatim}
yabby> seq_load cad.fas cad

 Reading the file 'cad.fas' ..
 6 sequence(s) found.

yabby> blast -E 5.0 tmp/LmjFwholegenome_20050601_V5.2.pep cad

 6 sequence(s) found in the object 'cad'

 Now running BLAST ..
 BLASTing sequence 1 of 6 (Q45153_BACFI)
 Query sequence 'Q45153_BACFI'
 Database 'tmp/LmjFwholegenome_20050601_V5.2.pep'
 Found 2 hits above the threshold (E=5.00)
 The best hit: 'LmjF36.2210'
 E-score = 1.93e+00

 BLASTing sequence 2 of 6 (Q53650_STAAU)
 Query sequence 'Q53650_STAAU'
 Database 'tmp/LmjFwholegenome_20050601_V5.2.pep'
 Found 1 hits above the threshold (E=5.00)
 The best hit: 'LmjF04.0510'
 E-score = 3.57e+00

 BLASTing sequence 3 of 6 (Q97PJ0_STRPN)
 Query sequence 'Q97PJ0_STRPN'
 Database 'tmp/LmjFwholegenome_20050601_V5.2.pep'
 Found 1 hits above the threshold (E=5.00)
 The best hit: 'LmjF13.1210'
 E-score = 2.11e+00

 BLASTing sequence 4 of 6 (P95773_STALU)
 Query sequence 'P95773_STALU'
 Database 'tmp/LmjFwholegenome_20050601_V5.2.pep'
 No BLAST hits above the threshold (E=5.00) found.

 BLASTing sequence 5 of 6 (Q9JXN6_NEIMB)
 Query sequence 'Q9JXN6_NEIMB'
 Database 'tmp/LmjFwholegenome_20050601_V5.2.pep'
 No BLAST hits above the threshold (E=5.00) found.

 BLASTing sequence 6 of 6 (Q9CE92_LACLA)
 Query sequence 'Q9CE92_LACLA'
 Database 'tmp/LmjFwholegenome_20050601_V5.2.pep'
 No BLAST hits above the threshold (E=5.00) found.
\end{verbatim}

The BLAST search was executed for each sequences in the 'cad' set
in turn, creating six 'blast' objects: 

\begin{verbatim}
yabby> what

    object(s)      type
  ------------------------------
    cad            seq           
    cad_1          blast         
    cad_2          blast         
    cad_3          blast         
    cad_4          blast         
    cad_5          blast         
    cad_6          blast         
\end{verbatim}

These can be examined later with the command 'blast\_info':

\begin{verbatim}
yabby> blast_info cad_1

 Query sequence 'Q45153_BACFI'
 Database 'tmp/LmjFwholegenome_20050601_V5.2.pep'
 Found 2 hits above the threshold (E=5.00)
 The best hit: 'LmjF36.2210'
 E-score = 1.93e+00
\end{verbatim}

\section{EMBOSS 'needle' commands}

\index{emboss}
\index{emboss\_needle}

The command 'emboss\_needle' shows the sequences from the output
of the program 'needle' from the EMBOSS suite of programs
\cite{emboss}. The input to the program 'needle' are two
sequences, and the program calculates the global sequence
alignment (Needleman-Wunsch algorithm). Typically two files
are provided as the input, each containing one sequence in FASTA  
format. If the second file contains multiple sequences,
the global sequence alignment will be calculate between
the first sequence and each sequence in turn from the
second set. If the second set contains many sequences, the
output file will be large, and it may be difficult to see
the best matches. These can be extracted with the command
'emboss\_needle':

\begin{verbatim}
yabby> emboss_needle needle.out 3

 Processing the file 'needle.out' ..

 (1) Q45153_BACFI:LmjF05.1170, Similarity: 26.6
 (2) Q45153_BACFI:LmjF05.1120, Similarity: 16.4
 (3) Q45153_BACFI:LmjF05.1040, Similarity: 15.0
\end{verbatim}

The second argument is the number of top-scoring sequences
to be printed out.

\index{emboss\_needl2}

The purpose of the command 'emboss\_needl2' is to find out
what is the similarity between the two sets of sequences.
Consider the following problem: given two sets of sequences,
set A and set B, we need to find out how similar are the
sequences in the set B to the sequences in the set A.

The command emboss\_needl2 will take one sequence at a time
from the query set A, and run the EMBOSS program 'needle' to
find the similarity after the best alignment between this
sequence and each sequence in the set B. emboss\_needl2 will
then extract the ID of the most similar sequence, as reported
by 'needle', and it will move to the next query sequence from
the query set A set until all query sequences are exausted.

Consider the following example: the set A is 'Ecoli057.fas'
and the set B is 'Ecoli\_lab.fas', both given as FASTA files.
To find out the similarity of sequences in Ecoli057.fas
to sequences in Ecoli\_lab.fas, we first load the query set
in yabby:

\begin{verbatim}
yabby> seq_load Ecoli057.fas Ecoli057

 Reading the file 'Ecoli057.fas' ..
 9 sequence(s) found.
\end{verbatim}

The second set 'Ecoli\_lab.fas' is required as FASTA file.
Then:

\begin{verbatim}
yabby> emboss_needl2 Ecoli057 Ecoli_lab.fas sim_list

 'Ecoli057' contains 9 sequence(s)
 Working on sp|P36546|OMPX_ECOLI
Needleman-Wunsch global alignment.
 Working on tr|Q8X8K6|Q8X8K6
Needleman-Wunsch global alignment.
 Working on sp|P39170|UP05_ECOLI
Needleman-Wunsch global alignment.
 Working on tr|Q8XDF1|Q8XDF1
Needleman-Wunsch global alignment.
 Working on sp|P58603|OMPT_ECO57
Needleman-Wunsch global alignment.
 Working on tr|Q8XAS0|Q8XAS0
Needleman-Wunsch global alignment.
 Working on tr|Q8X592|Q8X592
Needleman-Wunsch global alignment.
 Working on tr|Q8X7N5|Q8X7N5
Needleman-Wunsch global alignment.
 Working on tr|Q8X8F2|Q8X8F2
Needleman-Wunsch global alignment.
yabby>
\end{verbatim}

The argument to the command  'emboss\_needl2' are the query sequence
object, the FASTA file of the database set, and the name of the
object to be created which contains the information about similarities
This is 'needl2' object:
 
\begin{verbatim}
yabby> what

    object(s)      type
  ------------------------------
    Ecoli057       seq
    sim_list       needl2
\end{verbatim}

One can print the 'needl2' object:

\begin{verbatim}
yabby> print sim_list.needl2

OMPX_ECOLI OMPX_ECOLI 100.0
Q8X8K6 OMPG_ECOLI 100.0
UP05_ECOLI UP05_ECOLI 100.0
Q8XDF1 OMPF_ECOLI 99.4
OMPT_ECO57 OMPT_ECOLI 99.7
Q8XAS0 NMPC_ECOLI 80.3
Q8X592 MIPA_ECOLI 99.6
Q8X7N5 PHOE_ECOLI 99.7
Q8X8F2 YFEN_ECOLI 32.1
\end{verbatim}

Or sorted by percent similarity:

\begin{verbatim}
yabby> print -s sim_list.needl2

Q8X8F2 YFEN_ECOLI 32.1
Q8XAS0 NMPC_ECOLI 80.3
Q8XDF1 OMPF_ECOLI 99.4
Q8X592 MIPA_ECOLI 99.6
OMPT_ECO57 OMPT_ECOLI 99.7
Q8X7N5 PHOE_ECOLI 99.7
OMPX_ECOLI OMPX_ECOLI 100.0
Q8X8K6 OMPG_ECOLI 100.0
UP05_ECOLI UP05_ECOLI 100.0
\end{verbatim}

The above table lists two sequence IDs (ID\_query and ID\_dba) and
the percent similarity as reported by the program 'needle' after
the best global alignment.  The above table contains nine rows,
since the query set ('Ecoli057.fas') contains nine sequences.  Each
line refers to one sequence from the 'Ecoli057.fas', and given are
the best matching sequence from the second set ('Ecoli\_lab.fas'),
and percent similarity. We see that the sequence 'Q8X8F2' from
'Ecoli057.fas' has the least similarity to any sequence from the
'Ecoli\_lab.fas' set.

The table produced by 'emboss\_needl2' can be written to a file:

\begin{verbatim}
yabby> print -f needl2.out -s sim_list.needl2

 'sim_list.needl2' written to the file 'needl2.out'
\end{verbatim}


\section{Stand-alone commands}

Stand-alone commands do not interact with other commands. They
are, in effect, independent program scripts.  Nevertheless it is
often useful to collect such independent scripts under Yabby
to simplify their maintenance and ensure long-term retention.
An example of such command is 'pdb\_model'.

\index{pdb\_model.pl}

Another example of a stand-alone command is 'pdb\_model', which 
splits the Protein Data Bank (PDB) file with multiple models.
The output is a series of PDB files, one per model:

\begin{verbatim}
yabby> pdb_model 2C7H.pdb RBP_

 working on model  1 (creating 'RBP_1.pdb')
 working on model  2 (creating 'RBP_2.pdb')
 working on model  3 (creating 'RBP_3.pdb')
 working on model  4 (creating 'RBP_4.pdb')
 working on model  5 (creating 'RBP_5.pdb')
.... extra output deleted ....
\end{verbatim}

where '2C7H.pdb' is the original PDB file with multiple models,
and 'RBP\_' is the prefix used for individual model files.

\section{Extending Yabby}

The one-to-one correspondence between Yabby commands and Yabby
command scripts makes it very clear from where the particular
funcionality comes from. This also encapsulates a functionality
within a particular script file, and small size of individual
script files makes it easy to understand, modify, and maintain
such scripts.  The complete decoupling between command scripts
ensures that it is hard to break existing commands when adding
new commands.

Any functional Perl script placed in the Yabby library directory
will become available for execution as Yabby command.  Consider
the file 'hello.pl' with the following content:

\begin{verbatim}
print "Hello World!\n";
\end{verbatim}

If placed in the the Yabby library directory, it becomes immediately
available as the command 'hello':

\begin{verbatim}
$ yabby

 - YABBY version 0.1 - 

 [ 36 command(s) ready ]

yabby> hello

Hello World!

yabby>
\end{verbatim}

If already running, Yabby will need to be restarted before the command
'hello' becomes available.

\subsection{Creating a new Yabby command: seq\_letter}

Although a simple perl script would work if just dropped in the Yabby
library, to create a Yabby command would require some additional effort.
Consider creating a very simple command which operates on the sequence
object, and prints the first letter of the first sequence stored in
the sequence object. The command will be named 'seq\_letter'.

We start by creating the file 'seq\_letter.pl' in the Yabby library:

\begin{verbatim}
# seq_letter.pl

use yabby_sys;
use yabby_seq;

use Getopt::Std;
\end{verbatim}

All Yabby commands rely on the module 'yabby\_sys'. The module 'yabby\_seq'
contains functions for the manipulation of sequences. The Perl module
'Getopt::Std' is required to enable processing of single-character
command switches.

The next step is to create the usage string:

\begin{verbatim}
$USAGE = "
 Prints the first letter from the first sequence of a sequence
 object. 

 Usage:
        seq_letter OBJ_NAME

 where OBJ_NAME is the name of an existing sequence object.

 Notes:

 1. This command is merely an example to demonstrate how new
 Yabby commands can be developed.
";
\end{verbatim}

The usage string will be printed when 'help seq\_letter' or
'seq\_letter help' is issued.

We add the initialization commands:

\begin{verbatim}
# options
# initialization
@argl = sys_init( @ARGV );
check_call( @argl, [ 1 ] );
$obj_name = $argl[0];
\end{verbatim}

In the case of this command there are no options. The function
'sys\_init' is used to process arguments to all Yabby commands
in a standard way. The function 'check\_call' is used to check
the number of arguments.  In this case, exactly one argument
must be given, the name of an existing sequence object.  Finally,
in the last line the name of sequence object is associated to
the variable 'obj\_name'. In the next step the sequence object
is loaded in the memory:
 
\begin{verbatim}
requirements( $obj_name, $SEQUENCE );
$xmldoc = load_ip_xml( $obj_name, $SEQUENCE );
\end{verbatim}

The function 'requirements()' checks whether the sequence object
named 'obj\_name' exists. If not, it will terminate the script
with an appropriate error message. If this step is passed, the
function 'load\_ip\_xml()' loads the sequence object from the
workspace. Finally, we fetch the first sequence and print its
first letter:

\begin{verbatim}
$seq_hash = xml2seq( $xmldoc );
$seq_item = $seq_hash->{1};
printf " The first letter of the first sequence is: '%s'\n",
  substr( $seq_item->{$DBA_SEQUENCE}, 0, 1 );
\end{verbatim}

The complete listing of the modified command script 'seq\_letter'
is given below: 

\begin{verbatim}
# seq_letter.pl

use yabby_sys;
use yabby_seq;

use Getopt::Std;

$USAGE = "
 Prints the first letter from the first sequence of a sequence
 object. 

 Usage:
        seq_letter OBJ_NAME

 where OBJ_NAME is the name of an existing sequence object.

 Notes:

 1. This command is merely an example to demonstrate how new
 Yabby commands can be developed.
";

# options
# initialization
@argl = sys_init( @ARGV );
check_call( @argl, [ 1 ] );
$obj_name = $argl[0];

# requirements
requirements( $obj_name, $SEQUENCE );
$xmldoc = load_ip_xml( $obj_name, $SEQUENCE );

# body
$seq_hash = xml2seq( $xmldoc );
$seq_item = $seq_hash->{1};
printf " The first letter of the first sequence is: '%s'\n",
  substr( $seq_item->{$DBA_SEQUENCE}, 0, 1 );
\end{verbatim}

When executed, the command gives the following output:

\begin{verbatim}
yabby> seq_load tom20.fas tom20

 Reading the file 'tom20.fas' ..
 3 sequence(s) found.

yabby> seq_letter tom20

 The first letter of the first sequence is: 'M'
\end{verbatim}

Consider what happens when no arguments are given:

\begin{verbatim}
yabby> seq_letter 

 seq_letter:: ERROR: exactly 1 argument(s) required [ 0 supplied ]
 [ command 'seq_letter' failed ]
\end{verbatim}

This error was raised by the function 'check\_call'. If the sequence
object does not exist:

\begin{verbatim}
yabby> seq_letter ttt

 seq_letter:: missing requirement: 'ttt.seq'

 [ command 'seq_letter' failed ]
\end{verbatim}

The last error was raised by the function 'requirements()'.

\subsection{Extending seq\_letter to work on all sequences}

It would be easy to extend the command seq\_letter to print the
first letter of all sequences from the sequence object. This would
require modifying the body of the command:

\begin{verbatim}
# body
$seq_hash = xml2seq( $xmldoc );
$keys = get_seq_keys( $seq_hash );
printf " '%s' contains %d sequence(s)\n", $obj_name, $#{$keys}+1;

for $key ( @$keys ) {
  $seq_item = $seq_hash->{$key};
  printf " Sequence '%s', the first letter is '%s'\n",
  $seq_item->{$DBA_SEQID}, substr( $seq_item->{$DBA_SEQUENCE}, 0, 1 );
}
\end{verbatim}

The output of 'seq\_letter' is:

\begin{verbatim}
yabby> seq_letter tom20

 'tom20' contains 3 sequence(s)
 Sequence 'A.thaliana2', the first letter: 'M'
 Sequence 'O.sativa', the first letter: 'M'
 Sequence 'L.esculentum', the first letter: 'M'
\end{verbatim}

\subsection{Adding a command switch}

Consider how one could add an optional switch to print the last
letter of each sequence, instead of the first letter.

To add the switch '-t':

\begin{verbatim}
# options
getopts('t');

if ( defined( $opt_t ) ) {
  $opt_t_flag = 1;
} else {
  $opt_t_flag = 0;
}
\end{verbatim}

To extend the command body to take into account possible '-t'
switch: 

\begin{verbatim}
$seq_hash = xml2seq( $xmldoc );
$keys = get_seq_keys( $seq_hash );
printf " '%s' contains %d sequence(s)\n", $obj_name, $#{$keys}+1;

for $key ( @$keys ) {
  $seq_item = $seq_hash->{$key};
  $seq_id = $seq_item->{$DBA_SEQID};
  $sequence = $seq_item->{$DBA_SEQUENCE};

  if ( $opt_t_flag ) {
    $lett = substr( $seq_item->{$DBA_SEQUENCE}, -1, 1 );
  } else {
    $lett = substr( $seq_item->{$DBA_SEQUENCE}, 0, 1 );
  }

  printf " Sequence '%s', the letter is '%s'\n", $seq_id, $lett;
}
\end{verbatim}

The possible use of the modified 'seq\_letter':

\begin{verbatim}
yabby> seq_letter tom20

 'tom20' contains 3 sequence(s)
 Sequence 'A.thaliana2', the letter is 'M'
 Sequence 'O.sativa', the letter is 'M'
 Sequence 'L.esculentum', the letter is 'M'

yabby> seq_letter -t tom20

 'tom20' contains 3 sequence(s)
 Sequence 'A.thaliana2', the letter is 'R'
 Sequence 'O.sativa', the letter is 'R'
 Sequence 'L.esculentum', the letter is 'R'
\end{verbatim}

The full listing of the command 'seq\_letter' is:

\begin{verbatim}
# seq_letter.pl

use yabby_sys;
use yabby_seq;

use Getopt::Std;

$USAGE = "
 Prints the first sequence letter from the sequence object. 

 Usage:
        seq_letter [ options ] OBJ_NAME

 where OBJ_NAME is the name of an existing sequence object.

 Options:

 -t -- print the last letter instead of the first

 Notes:

 1. This command is merely an example to demonstrate how new
 Yabby commands can be developed.
";

# options
getopts('t');

if ( defined( $opt_t ) ) {
  $opt_t_flag = 1;
} else {
  $opt_t_flag = 0;
}

# initialization
@argl = sys_init( @ARGV );
check_call( @argl, [ 1 ] );
$obj_name = $argl[0];

# requirements
requirements( $obj_name, $SEQUENCE );
$xmldoc = load_ip_xml( $obj_name, $SEQUENCE );

# body
$seq_hash = xml2seq( $xmldoc );
$keys = get_seq_keys( $seq_hash );

printf " '%s' contains %d sequence(s)\n", $obj_name, $#{$keys}+1;

for $key ( @$keys ) {

  $seq_item = $seq_hash->{$key};
  $seq_id = $seq_item->{$DBA_SEQID};
  $sequence = $seq_item->{$DBA_SEQUENCE};

  if ( $opt_t_flag ) {
    $lett = substr( $seq_item->{$DBA_SEQUENCE}, -1, 1 );
  } else {
    $lett = substr( $seq_item->{$DBA_SEQUENCE}, 0, 1 );
  }

  printf " Sequence '%s', the letter is '%s'\n", $seq_id, $lett;
}
\end{verbatim}

\subsection{A command that creates an object}

The command 'seq\_letter' is very simple in that it merely prints
something on the output screen; it does not modify an existing
object nor does it creates any new Yabby objects. Consider a command
that actually creates a sequence object. An example for this may
be a command that chops the first letter from each sequence in
the sequence object, and saves the result as a new sequence object.
This requires two names to be given on the input, the name of
an existing sequence object, and the name of the sequence object
that will be created:

\begin{verbatim}
seq_letterc OBJ_NAME OBJ_NAME_NEW
\end{verbatim}

The new command is named 'seq\_letterc' to distinguish it from
the previous 'seq\_letter'.

To build the command 'seq\_letterc' we create the file 'seq\_letterc.pl'
in the Yabby library:

\begin{verbatim}
# seq_letters.pl

use yabby_sys;
use yabby_seq;

use Getopt::Std;

$USAGE = "
 Chops the first sequence letter from the sequence object, and
 saves the result as a new sequence object.

 Usage:
        seq_letters OBJ_NAME OBJ_NAME_NEW

 where OBJ_NAME is the name of an existing sequence object.

 Notes:

 1. This command is merely an example to demonstrate how new
 Yabby commands can be developed.
";
\end{verbatim}

Two arguments are now required:

\begin{verbatim}
# options
# initialization
@argl = sys_init( @ARGV );
check_call( @argl, [ 2 ] );
$obj_name = $argl[0];
$obj_name2 = $argl[1];
\end{verbatim}

We can use the same construct to check for requirements and fetch the
sequence object:

\begin{verbatim}
# requirements
requirements( $obj_name, $SEQUENCE );
$xmldoc = load_ip_xml( $obj_name, $SEQUENCE );

# body
$seq_hash = xml2seq( $xmldoc );
$keys = get_seq_keys( $seq_hash );

printf " '%s' contains %d sequence(s)\n", $obj_name, $#{$keys}+1;
\end{verbatim}

The new command requires a new sequence object to be created, and
populated with sequences from the existing sequence object (minus
the first letter). We first initialize a new sequence object:

\begin{verbatim}
$seq_hash2 = {};
\end{verbatim}

The same construct as in 'seq\_letter' could be used to loop through
the sequences. Now, we also create a new sequence for each existing
sequence, chop the first letter, and joing the result to the new
sequence object:

\begin{verbatim}
for $key ( @$keys ) {

  $seq_item = $seq_hash->{$key};

  $seq_item_new = {};
  $seq_item_new->{$DBA_SEQID} = $seq_item->{$DBA_SEQID};
  $seq_item_new->{$DBA_COMMENT} = $seq_item->{$DBA_COMMENT};
  $seq_item_new->{$DBA_SEQUENCE} = $seq_item->{$DBA_SEQUENCE};

  # chop the first letter from the sequence
  substr( $seq_item_new->{$DBA_SEQUENCE}, 0, 1 ) = "";

  $seq_hash2->{$key} = $seq_item_new;
}
\end{verbatim}

We need to save the newly created sequence hash 'seq\_hash2' as the
new sequence object. In the first step, we convert the sequence hash
object into an XML document:

\begin{verbatim}
$xmldoc = seq2xml( $seq_hash2 );
\end{verbatim}

This is the reverse of the above command 'xml2seq()', which converts
an XML document representing a sequence object into a sequence hash,
the form all commands employ to represent the sequence object in
memory.  The functions 'xml2seq()' and 'seq2xml()' are specific for
sequence objects, and are defined in the package 'yabby\_seq.pm'.
Finally, we store the XML sequence object in the workspace:

\begin{verbatim}
printf " Saving chopped sequences as '%s'\n", $obj_name2;
save_ip_xml( $xmldoc, $obj_name2, $SEQUENCE, $WARN_OVERW );
\end{verbatim}

The function 'save\_ip\_xml()' sends the XML representation of an
object to the workspace. This function is the reverse analog of
the previously used function 'load\_ip\_xml()', which fetches the
XML representation of an object from the workspace. These two
functions are defined in the package 'yabby\_sys.pm'.

The arguments to the function 'save\_ip\_xml()' are are follows:
the XML representation of an object ('xmldoc'), the name of the
object to be created ('obj\_name2'), the type of the object to be
created ('SEQUENCE' in this case), and the flag which determines
whether the warning will be issued is the object with the same
name and the same type already exists in the workspace and will
be overwritten ('WARN\_OVERW').

The command 'seq\_letterc' could be used as follows:

\begin{verbatim}
yabby> seq_load tom20.fas tom20

 Reading the file 'tom20.fas' ..
 3 sequence(s) found.

yabby> what      

    object(s)      type
  ------------------------------
    tom20          seq           

yabby> seq_letterc tom20 tom20_chop

 'tom20' contains 3 sequence(s)
 Saving chopped sequences as 'tom20_chop'

yabby> what

    object(s)      type
  ------------------------------
    tom20          seq           
    tom20_chop     seq           

yabby> print tom20.seq

>A.thaliana2 [ A.thaliana2 ]
MEFSTADFERFIMFEHARKNSEAQYKNDPLDSENLLKWGGALLELSQFQPIPEAKLMLND
AISKLEEALTINPGKHQALWCIANAYTAHAFYVHDPEEAKEHFDKATEYFQRAENEDPGN
DTYRKSLDSSLKAPELHMQFMNQGMGQQILGGGGGGGGGGMASSNVSQSSKKKKRNTEFT
YDVCGWIILACGIVAWVGMAKSLGPPPPAR
>O.sativa [ O.sativa ]
MDMGAMSDPERMFFFDLACQNAKVTYEQNPHDADNLARWGGALLELSQMRNGPESLKCLE
DAESKLEEALKIDPMKADALWCLGNAQTSHGFFTSDTVKANEFFEKATQCFQKAVDVEPA
NDLYRKSLDLSSKAPELHMEIHRQMASQASQAASSTSNTRQSRKKKKDSDFWYDVFGWVV
LGVGMVVWVGLAKSNAPPQAPR
>L.esculentum [ L.esculentum ]
MDMQSDFDRLLFFEHARKTAETTYATDPLDAENLTRWAGALLELSQFQSVSESKKMISDA
ISKLEEALEVNPQKHDAIWCLGNAYTSHGFLNPDEDEAKIFFDKAAQCFQQAVDADPENE
LYQKSFEVSSKTSELHAQIHKQGPLQQAMGPGPSTTTSSTKGAKKKSSDLKYDVFGWVIL
AVGLVAWIGFAKSNMPXPAHPLPR

yabby> print tom20_chop.seq

>A.thaliana2 [ A.thaliana2 ]
EFSTADFERFIMFEHARKNSEAQYKNDPLDSENLLKWGGALLELSQFQPIPEAKLMLNDA
ISKLEEALTINPGKHQALWCIANAYTAHAFYVHDPEEAKEHFDKATEYFQRAENEDPGND
TYRKSLDSSLKAPELHMQFMNQGMGQQILGGGGGGGGGGMASSNVSQSSKKKKRNTEFTY
DVCGWIILACGIVAWVGMAKSLGPPPPAR
>O.sativa [ O.sativa ]
DMGAMSDPERMFFFDLACQNAKVTYEQNPHDADNLARWGGALLELSQMRNGPESLKCLED
AESKLEEALKIDPMKADALWCLGNAQTSHGFFTSDTVKANEFFEKATQCFQKAVDVEPAN
DLYRKSLDLSSKAPELHMEIHRQMASQASQAASSTSNTRQSRKKKKDSDFWYDVFGWVVL
GVGMVVWVGLAKSNAPPQAPR
>L.esculentum [ L.esculentum ]
DMQSDFDRLLFFEHARKTAETTYATDPLDAENLTRWAGALLELSQFQSVSESKKMISDAI
SKLEEALEVNPQKHDAIWCLGNAYTSHGFLNPDEDEAKIFFDKAAQCFQQAVDADPENEL
YQKSFEVSSKTSELHAQIHKQGPLQQAMGPGPSTTTSSTKGAKKKSSDLKYDVFGWVILA
VGLVAWIGFAKSNMPXPAHPLPR
\end{verbatim}

Since no safeguards were put against overwriting an existing
sequence object, the command 'seq\_letterc' could be used
twice to chop sequences from the N-terminus in place,

\begin{verbatim}
yabby> seq_letterc tom20 tom20

 'tom20' contains 3 sequence(s)
 Saving chopped sequences as 'tom20'
 [ seq_letterc: 'tom20.seq' exists, overwritten ]

yabby> seq_letterc tom20 tom20

 'tom20' contains 3 sequence(s)
 Saving chopped sequences as 'tom20'
 [ seq_letterc: 'tom20.seq' exists, overwritten ]
\end{verbatim}

\subsection{Creating a new object type}

Yabby supports two internal representations of data in the workspace:
simple lists and XML objects. To create a new object type one must
first decide whether the data in question can be effectively
represented with a two-dimensional table. Two-dimensional tables
can be effectively represented with two-dimensional lists, where
each item in the list corresponds to the table row. Lists are simple
and efficient in Perl, and if suitable this is the preferred
representation of data. However, XML format is much more powerful
and flexible, and is the preferred solution when the data has
a complex internal structure.

A sequence object uses XML to store data into the workspace.
A Yabby commands store sequence objects
in memory in the form of a hash (sequence hash).  All sequence
commands use the standard function 'load\_ip\_xml()' to fetch the XML
representation of a sequence object from the workspace. Subsequently
the object is converted into a more efficient structure, a
sequence hash, with the function 'xml2seq()'. After the sequences
are manipulated, or a new sequence object is created, the result
may be stored in the workspace. In this case, the reverse
transformation process takes place:  the sequence hash is convereted
into an XML object with the function 'seq2xml()', and the resulting
XML object is saved in the workspace with the function 'save\_ip\_xml()'.

The functions 'load\_ip\_xml()' and 'save\_ip\_xml()' can be used
to load/store the XML representation of any object in the workspace.
The analogous functions for objects represented as a two-dimensional
lists are 'load\_ip()' and 'save\_ip()'.

As an example we will create a new object type named "seqid" to
represent the list of sequence IDs. To demonstrate the use of this
object we will also create the command 'seq2id' which creates the
"seqid" object from an existing sequence object.

The first step is to create the script 'seq2id.pl' in the Yabby
library, and we start as previously:

\begin{verbatim}
# seq2id.pl

use yabby_sys;
use yabby_seq;

use Getopt::Std;

$USAGE = "
 Converts a sequence object into seqid object.

 Usage:
        seq2id OBJ_NAME

 where OBJ_NAME is the name of an existing sequence object.

 Notes:

 1. This command is merely an example to demonstrate how a
 new object type can be implemented.
";

# options
# initialization
@argl = sys_init( @ARGV );
check_call( @argl, [ 1 ] );
$obj_name = $argl[0];

# requirements
requirements( $obj_name, $SEQUENCE );
$xmldoc = load_ip_xml( $obj_name, $SEQUENCE );

# body
$seq_hash = xml2seq( $xmldoc );
$keys = get_seq_keys( $seq_hash );
\end{verbatim}

With this we have the sequence hash in memory. In 'yabby\_seq.pm' we
add the name of the "seqid" object:

\begin{verbatim}
$SEQUENCE = "seq";
$MOTIF = "motif";
$SEQID = "seqid";
\end{verbatim}

The new object will be represented as a one-dimensional list, where
each sequence ID string will be one element of this list. This is
equivalent to a two-dimensional table, where there is only one
column in the table. Next we create the variable 'seqid\_obj' to
hold the 'seqid' object, fill this object with the values from
the sequence hash, and store the object in the workspace by using
the standard function:

\begin{verbatim}
$seqid_obj = [];

for $key ( @$keys ) {

  $seq_item = $seq_hash->{$key};
  $seq_id = $seq_item->{$DBA_SEQID};

  push @$seqid_obj, [ $seq_id ];
}

print " Saving '$obj_name.$SEQID'\n";
save_ip( $seqid_obj, $obj_name, $SEQID, $WARN_OVERW );
\end{verbatim}

The complete listing for the command script 'seq2id':

\begin{verbatim}
# seq2id.pl

use yabby_sys;
use yabby_seq;

use Getopt::Std;

$USAGE = "
 Converts a sequence object into seqid object.

 Usage:
        seq2id OBJ_NAME

 where OBJ_NAME is the name of an existing sequence object.

 Notes:

 1. This command is merely an example to demonstrate how a
 new object type can be implemented.
";

# options
# initialization
@argl = sys_init( @ARGV );
check_call( @argl, [ 1 ] );
$obj_name = $argl[0];

# requirements
requirements( $obj_name, $SEQUENCE );
$xmldoc = load_ip_xml( $obj_name, $SEQUENCE );

# body
$seq_hash = xml2seq( $xmldoc );
$keys = get_seq_keys( $seq_hash );

printf " '%s' contains %d sequence(s)\n", $obj_name, $#{$keys}+1;

$seqid_obj = [];

for $key ( @$keys ) {

  $seq_item = $seq_hash->{$key};
  $seq_id = $seq_item->{$DBA_SEQID};

  push @$seqid_obj, [ $seq_id ];
}

print " Saving '$obj_name.$SEQID'\n";
save_ip( $seqid_obj, $obj_name, $SEQID, $WARN_OVERW );
\end{verbatim}

And below is the output when the command is executed on an example
sequence object:

\begin{verbatim}
yabby> seq_load tom20.fas tom20

 Reading the file 'tom20.fas' ..
 3 sequence(s) found.

yabby> seq2id tom20

 'tom20' contains 3 sequence(s)
 Saving 'tom20.seqid'

yabby> what

    object(s)      type
  ------------------------------
    tom20          seq           
                   seqid         
\end{verbatim}

Of course, at this point nothing can be done with the 'seqid' object.
It cannot be even printed on the screen:

\begin{verbatim}
yabby> print tom20.seqid

 print:: ERROR: printing the property 'seqid' not yet implemented
 [ command 'print' failed ]
\end{verbatim}

This is because the printing of each particular object type has to
be enabled in the 'print' command. However, one can view the 'seqid'
object directly from the filesystem, since 'tom20.seqid' is just
a simple text file, containing sequence IDs listed one per row:

\begin{verbatim}
yabby> cat .yabby/tom20.seqid
A.thaliana2
O.sativa
L.esculentum
\end{verbatim}

\subsection{Some additional explanations and programming conventions}

All objects that are represents with two dimensional tables are
initialized as annonymous arrays:

\begin{verbatim}
$seqid_obj = [];
\end{verbatim}

When filled with values, each element of this array is a reference
to anonymous array. In this case containing only one element 'seq\_id':

\begin{verbatim}
push @$seqid_obj, [ $seq_id ];
\end{verbatim}

This format is very flexible.  If the data table needs to contain more
than one element per row, these would be simply added to the anonymous
reference representing a table row:

\begin{verbatim}
push @$seqid_obj, [ $seq_id, $colum2_elem, $colum3_elem, ... ];
\end{verbatim}

Consider how the object 'mol' is built, which is internally represented
with a two-dimensional table:

\begin{verbatim}
push @$mol, [ $resi_num, $resi_name, $atom_name, $segmentID ];
\end{verbatim}

The function 'save\_ip()' expects the reference to an anonymous array,
where each element of this array is a reference to another anonymous
array which may contain one or more elements. This holds for any
object represented by a list, and functions 'save\_ip()' and
'load\_ip()' are used to store/fetch the object from the workspace.
For objects represented in XML, a special converter function needs
to be written akin to 'xml2seq()' and 'seq2xml()'.
